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Over the last decade, measurements of the CMB anisotropy has spearheaded the remark- 
able transition of cosmology into a precision science. However, addressing the systematic 
' effects in the increasingly sensitive, high resolution, 'full' sky measurements from different 

• CMB experiments pose a stiff challenge. The analysis techniques must not only be computa- 

tionally fast to contend with the huge size of the data, but, the higher sensitivity also limits 
Q^. the simplifying assumptions which can then be invoked to achieve the desired speed without 

compromising the final precision goals. While maximum likelihood is desirable, the enormous 
computational cost makes the suboptimal method of power spectrum estimation using Pseudo- 
^ ' C; unavoidable for high resolution data. The debiasing of the Pseudo-C; needs account for 

Ci ■ non-circular beams, together with non-uniform sky coverage. We provide a (semi)analytic 

framework to estimate bias in the power spectrum due to the effect of beam non-circularity 
and non-uniform sky coverage including incomplete/masked sky maps and scan strategy. The 
' approach is perturbative in the distortion of the beam from non-circularity, allowing for rapid 

^ . computations when the beam is mildly non-circular. We advocate that it is computationally 

advantageous to employ 'soft' azimuthally apodized masks whose spherical harmonic trans- 
form die down fast with m. We numerically implement our method for non-rotating beams. 
We present preliminary estimates of the computational cost to evaluate the bias for the up- 
coming CMB anisotropy probes (^max ^ 3000), with angular resolution comparable to the 
Planck surveyor mission. We further show that this implementation and estimate is applica- 
ble for rotating beams on equal declination scans and possibly can be extended to simple 
approximations to other scan strategies. 
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1 INTRODUCTION 

The fluctuations in the Cosmic Microwave Background (CMB) radiation are theoretically very well understood, allowing precise and un- 
ambiguous predictions for a given cosmological model (Bond 1996, Hu and Dodelson 2002). The measurement of CMB anisotropy with 
the ongoing Wilkinson Microwave Anisotropy Probe (WMAP) and the upcoming Planck surveyor, has ushered in a new era of precision 
cosmology. Such data rich experiments, with increased sensitivity, high resolution and 'full' sky measurements pose a stiff challenge for 
current analysis techniques to realize the full potential of precise determination of cosmological parameters. The analysis techniques must 
not only be computationally fast to contend with the huge size of the data, but, the higher sensitivity also limits the simplifying assumptions 
that can then be invoked to achieve the desired speed without compromising the final accuracy. As experiments improve in sensitivity, the 
inadequacy in modeling the observational reality start to limit the returns from these experiments. The current effort is to push the boundary 
of this inherent compromise faced by the current CMB experiments that measure the anisotropy in the CMB temperature and polarization. 
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Accurate estimation of tiie angular power spectrum, Ci, is inarguably tlie foremost concern of most CMB experiments. The extensive 
literature on this topic has been summarized (Hu and Dodelson 2002, Bond 1996, Efstathiou 2004). For Gaussian, statistically isotropic 
CMB sky, the Ci that corresponds to the covariance that maximizes the multivariate Gaussian PDF of the temperature map, Ar(q), is the 
Maximum Likelihood (ML) solution. Different ML estimators have been proposed and implemented on CMB data of small and modest 
sizes (Gorski 1994, Gorski et al. 1994, 1996, 1997, Tegmark 1997, Bond et al. 1998). While it is desirable to use optimal estimators of 
Ci that obtain (or iterate toward) the ML solution for the given data, these methods are usually limited by the computational expense of 
matrix inversion that scales as N'^ with data size Nd (Borrill 1999, Bond et al. 1999). Various strategies for speeding up ML estimation 
have been proposed, such as, exploiting the symmetries of the scan strategy (Oh et al. 1999), using hierarchical decomposition (Dore et al. 
2001), iterative multi-grid method (Pen 2003), etc. Variants employing linear combinations of AT(q) such as aj^ on set of rings in the sky 
can alleviate the computational demands in special cases (Challinor et al. 2002, van Leeuwen et al. 2002, Wandelt & Hansen 2003). Other 
promising 'exact' power estimation methods have been recently proposed(Knox et al. 2001, Jewell et al. 2004, Wandelt 2004). 

However there also exist computationally rapid, sub-optimal estimators of Ci. Exploiting the fast spherical harmonic transform (~ 
N^^^), it is possible to estimate the angular power spectrum Ci — \aim\'^ + 1) rapidly (Yu & Peebles 1969, Peebles 1974). This 
is commonly referred to as the Pseudo-C; method (Wandelt et al. 2003). Analogous approach employing fast estimation of the correlation 
function C(q • q') = {AT{q)AT{q')) have also been explored (Szapudi et al. 2001, Szapudi et al. 2001). It has been recently argued that 
the need for optimal estimators may have been over-emphasized since they are computationally prohibitive at large I. Sub-optimal estimators 
are computationally tractable and tend to be nearly optimal in the relevant high I regime. Moreover, already the data size of the current 
sensitive, high resolution, 'full sky' CMB experiments, such as WMAP have been compelled to use sub-optimal Pseudo-C; and related 
methods (Bennett et al. 2003, Hinshaw et al. 2003). On the other hand, optimal ML estimators can readily incorporate and account for 
various systematic effects, such as noise correlations, non-uniform sky coverage and beam asymmetries. A hybrid approach of using ML 
estimation of Ci for low I for low resolution map and Pseudo-C; like estimation for large / where it is nearly optimal has been suggested 
(Efstathiou 2004) and has even been employed by the recent analysis of the WMAP-3 year data (Hinshaw et al. 2007). The systematic 
correction to the Pseudo-C; power spectrum estimate arising from non-uniform sky coverage has been studied and implemented for CMB 
temperature (Hivon et al. 2002) and polarization (Brown at al 2005). The leading order systematic bias due to noncircular beam has been 
studied by us in an earlier publication (Mitra et al. 2004). Here we extend the results in a thorough analytical approach to include all the 
significant perturbation orders and combine the effect of non-uniform sky coverage. 

It has been usual in CMB data analysis to assume the experimental beam response to be circularly symmetric around the pointing 
direction. However, real beam response functions have deviations from circular symmetry. Even the main lobe of the beam response of 
experiments are generically non-circular (non-axisymmetric) since detectors have to be placed off-axis on the focal plane. (Side lobes and 
stray light contamination add to the breakdown of this assumption). For highly sensitive experiments, the systematic errors arising from the 
beam non-circularity become progressively more important. Dropping the circular beam assumption leads to major complications at every 
stage of the analysis pipeline. The extent to which the non-circularity of the beam affects the step of going from the time-stream data to sky 
map is very sensitive to the scan-strategy. The beam now has an orientation with respect to the scan path that can potentially vary along the 
path. This introduces an additional complication - that the beam function is inherently time dependent and difficult to deconvolve. 

Even after a sky map is made, the non-circularity of the effective beam affects the estimation of the angular power spectrum, Cj, by 
coupling the power at different multipoles, typically, on scales beyond the inverse angular beam-width. Mild deviations from circularity can 
be addressed by a perturbation approach (Souradeep & Ratra 2001, Fosalba et al. 2002) and the effect of non-circularity on the estimation 
of CMB power spectrum can be studied (semi) analytically (Mitra et al. 2004). Figure[T]shows the predicted level of bias due to noncircular 
beams in our formalism for elliptical beams compared to the noncircular beam corrections computed in the recent data release by WMAP 
(Hinshaw et al. 2007). 

To avoid contamination of the primordial CMB signal by Galactic emission, the region around the Galactic plane is masked from maps. 
If the Galactic cut is small enough, then the coupling matrix will be invertible, and the two-point correlation function can be determined 
on all angular scales from the data within the uncut sky (Mortlock et al. 2002). Hivon et al. (2002) present a technique (MASTER) for 
fast computation of the power spectrum accounting for the galactic cut, but restricted to circular beams. In our present work, we present 
analytical expressions for the bias matrix of the Pseudo-C; estimator for the incomplete sky coverage, using a noncircular beam. In Section|2] 
we show a heuristic approach to the analytic form of the bias matrix taking into account the above mentioned effects. We have shown that 
our estimation matches with the existing results in different limits in Section[3] In Section|4]we outline the numerical implementation of our 
approach, estimate the computational cost and suggest a potential algorithmic route to reducing the cost. The discussion and conclusion of 
this work is given in Section^ 



2 FORMALISM 

The CMB temperature anisotropy field Ar(q) over all the sky directions q = (6, 0) is assumed to be Gaussian and statistically isotropic 
and hence the angular power spectrum stores all the statistical information about the anisotropy field. These temperature fluctuations are 
expanded in the basis of spherical harmonics. 
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Figure 1. The predicted non-circular beam correction for a CMB experiment witli elliptical Gaussian beam with/iv/im beam-width of 0.51° and eccentricity 
e := \J\ — \P- la^ , where a (6) is the semi-major (semi-minor) axis, e = 0.65 for the Q beam (dashed line), beam-width of 0.35° and e = 0.46 for the 
V beam (dash-dotted line) and beam-width of 0.22° and e = 0.4 for the W beam (dash-dot-dot line) are shown in the left plot. The solid curves are the non 
circular beam corrections estimated by the WMAP team for the Q,V and W channels. In the plot in the right, we show that the non-circular beam correction 
matches with the estimates of the WMAP team for a slightly different change in the eccentricity of the Q and V beams (W beams matched in the left plot), 
with the new eccentricities e = 0.50 and e = 0.40 respectively. This difference is attributed to the fact that the beams visit the pixels multiple times with 
different orientations, and hence the effective eccentricity is reduced. 



AT(q) = ^a;„y,„(q), (1) 
where aim are the spherical harmonic transforms of the temperature anisotropy field 

aim ■- [ dne^AT{ci)Yim{ci). (2) 

The observed CMB temperature fluctuation field Ar(q) is a convolution of the "beam" profile i3(q, q') with the real temperature 
fluctuation field AT(q') and contaminated by additive experimental noise n(q). Moreover, due to the strong influence of the foreground 
emission in our galactic plane and point sources, some of the pixels have to be masked prior to power spectrum estimation. The mask function 
U (q) is usually assigned zero for the corrupt pixels and one for the clean ones, but it could be a smoother weight function also, that can take 
values between zero and one, as long as it sufficiently masks the foreground contamination. Mathematically the observed temperature can be 
expressed as 



AT(q) = [/(q) [/ dnq,B(q,q)Ar(q') + 7i(q) 
The spherical harmonic transform of the mask function. 



(3) 



Ulm = / dn^Yim{ei)U{ci), (4) 
Jin 

is a very useful quantity for our analysis. 

Statistical isotropy of underlying CMB anisotropy signal implies that the two point correlation function (7(qi, q2) = {AT(qi) Ar(q2)) 
depends only on the angular separation of the direction, i.e., C(qi, q2) ~ C(qi • q2). We can therefore expand it as a Fourier-Legendre 
series 

(AT(qi)AT(q2)) = J2 ^l^C7,P,(qi ■ q2), (5) 
where the coefficients C'l constitute the CMB power spectrum 

Ci = f d^e^, f dnq,(AT(qi)AT(q2))Pi(qi-q2)- (6) 

J 47r J 47r 

The addition theorem for spherical harmonics 

. I 

J2 YiU(li)Yim{fl2) = Pi(qi-q2), (7) 

m— — Z 

and the orthogonality of Legendre polynomials 

l'^dxPi{x)Pi,{x) = 2^5n' (8) 
can then be used to show that the matrix {aimaiim'} is diagonal: 
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{a-lmaVm') — ClSll'Smm'- (9) 

The observed two point correlation function for a statistically isotropic CMB anisotropy signal is 

C(q, q') = (AT (q)Ar(q')) = £ Ci Wi{c^, q') , (10) 

i=0 

where C; is the angular power spectrum of CMB anisotropy signal and the window function, 

H^i(qi,q2) := j dfiq j dt^q'B(qi,q)B(q2,q')J'i(q-q'). dD 

encodes the effect of finite resolution through the beam function. A CMB anisotropy experiment probes a range of angular scales characterized 
by this window function VK;(q, q'). The window depends both on the scanning strategy as well as the angular resolution and response of 
the experiment. However, it is neater to logically separate these two effects by expressing the window Wi{c{, q') as a sum of 'elementary' 
window functions of the CMB anisotropy (Souradeep & Ratra 2001). For a given scanning strategy, the results can be readily generalized 
using the representation of the window function as a sum over elementary window functions (see, e.g., Souradeep & Ratra 2001, Coble et al. 
2003, Mukherjee et al. 2003). 

For some experiments, the beam function may be assumed to be circularly symmetric about the pointing direction, i.e., _B(q, q') = 
_B(q ■ q'), without significantly affecting the results of the analysis. In any case, this assumption allows a great simplification, since the beam 
function can then be represented by an expansion in Legendre polynomials as 

B(q • q ) = ^ E (2« + 1) Bi p,(q ■ q ), (12) 

where, 

^ ■ ■ ■ (13) 



Bi = |\(q-q')Pi(q-q)e(q-q'), 



and Z?(q • q') is the circularized beam obtained by averaging B{z, q) over the azimuthal coordinate (f>. Consequently, it is straightforward to 
derive the well known simple expression 

Wi{ci, q') = Sf P(q-q'), (14) 
for a circularly symmetric beam function. For non-circular beams and incomplete sky coverage the Pseudo-Ci estimator, 

1 ' 

= 2iTi ^ ^^^^ 

m— — i 

takes the form 

Ci = ^ f dfle,, f dfiq, [/(qi)C/(q2)P(qi-q2)Af(qi)Af(q2). (16) 
In this case, the expectation value of the Pseudo-Ci estimator becomes 

(C,) = ^Aa.Cr + Cf, (17) 
I' 

that is, the estimator is non-trivially biased, where the bias matrix takes the form 

11'^ r ^ T 2 



^"' = ?7XtEE f dn^U{q)Yi„{(l) /dr!q,y;^(q')i3(q,q') 



(18) 



The noise term Cj^ , arising from instrumental noise, can be measured to a very high accuracy. If the noise term for full sky coverage Cj^ is 
known, it can be combined with the bias matrix for incomplete sky coverage A/;;/ (Hivon et al. 2002) to obtain the noise power spectrum for 
cut-sky 

=Y.MwC?, (19) 

V 

where, Af;;' is defined as 

with Wi" = \ Ui"m"\^ + !)■ Computation of the bias matrix is important for defining the unbiased Pseudo-C; estimator 
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(21) 



that removes the systematic effects of beam non-circularity and incomplete sky coverage. 

The experimental beams in CMB experiments are mildly non-circular, and hence it makes sense to define Beam Distortion Parameters 
(BDP) j3im{<t. 1) so that we can calculate the result in a perturbative expansion of small parameter. The BDP is expressed as = 
bim/bio, where the are the spherical harmonic moments of the beam function for the pointing direction z: 



d^nYil^iq) -B(z,q) , 

and Bi, defined in Eq. l ll3t . is connected to bio through the relation: 

r2TT 



Bi ^ I smf 

'0 



4-K 
21 + 1 



2tv 



.B(z,q) 



47r 



21 + 1 



(22) 



(23) 



Evaluation of the spherical harmonic transforms of the beam function -B(q, q') for each pointing direction q is computationally prohibitive. 
We use the spherical harmonic transforms of B{z, q'), incorporating rotation in it via Wigner-D functions, in order to compute the 
harmonic transforms of B(q, q') for any q by using the formula ( 142b 



/ dnq,y;,„(q')B(q,q) 

JiTT 



21' + 1 



^ Bi, D^^^, (q, /5(q)). 



m' — — I' 

Then, using the spherical harmonic expansion of the mask function U (q) 

oo I 

f/(q) = E E Uir^Yir^i.^)^ 

we can rewrite the general form of the bias matrix [Eq. dlSb l as 



Aw = 



where 



Bf {21' + 1) 
47r (2/ + 1) 



I i' 

E E 

■n— — l m — — l' 



oo I 



m" mm' 



l"^0 mf'^ — V' 



^nm" mm' 



mm. 



(24) 



(25) 



(26) 



(27) 



i+i" 



(-1) 



^(2; + l)(2Z" + l) 



47r 



L=\l-l"\ 



L+l' 



E^L (m — n — m ) /^L m L r /^N] 

^ L(-n-m")l'm ^ LOV -m' X{m-n-m" )m' [P(qj] 

L' = \L-l'\ 



and 



'[p(q)] ■= / cif^q-C'Lm'(q,p(q))- 

Ji-R 



(28) 



(29) 



To obtain the expression for the J coefficients above we have introduced Clebsch Gordon coefficients C^-^_^^^^, 
expansion of Wigner-d (see Appendix lAll for details). 

Eq. i28t provides the expression for the bias [Eq. | |26H in its full generality. For a specified general form of the scan strategy, however, 
we need to precompute the coefficients xlnm' which are functionals of p(q). Special cases, which often provide good approximations to the 
real scan strategies, provide significant computational advantages. From the definition of x coefficients we show below that 

I. For a /5(q) separable in declination and Right Ascension parts, i.e. p(q) = 0(6*) + 

xLm'[pm = r d^e-^-^^e-^""'"^^^ r sine dedl^, {9) e—''^^''^ (30) 
Jo Jo 

would have significant values only in a much constrained domain of the indices m,m' . 

II. For equal declination scan strategies, where p(q) = p(^'), we show below that the computational cost reduces to that of the non-rotating 
beam (p(q) — 0) case. 

Hence, the study of the bias computation for non-rotating beams provides an analytic framework that is computationally equivalent to a 
broader class of scan strategies as mentioned at the appropriate places below. 

It is clear from the definition of xlnm' l^hat for equal declination scans [case II] the (f) integral is leading to the constraint 



along with sinusoidal 
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Xmm' 



q) = 0] (X 

^0 



im(f> 



211 5mO- 



(31) 



So the J coefficients in Eq. ([28} contribute to bias only when m" — m — n, saving huge amount of computation. For a much broader class 
within case I, it is reasonable to expect that the x symbols in the final expression for bias will contribute when m" is close to m n. 
For equal declination scan strategies, the expression for the bias matrix can be written as 



Aw 



16tv 



•n=~-l m,— — l' 



J2 V2l" + 1 U, 



I" (m — n) 



l + l" 

'-■I0l"0 '^l 

L=\l-l"\ 



L + l' 



I' 



Lm ^L'O jL' 

lnl"{m-n) / , '^L-ml'm / ^ "OA 

L' = \L-l'\ N=-L' 

2 



L' [T^ 

■ON 



in' — — L' 

where the coefficients 

r„.'iv[p(e)] r' (-1)' 



n in iN9 —im p{6) 



(32) 



(33) 



are functionals of the scan strategy and have to be precomputed (analytically/numerically) for a given p(Q\ 

For an efficient computation of the J coefficients using a numerical implementation scheme described in section |4l we derive an 
alternative expression for the J coefficients using the sinusoidal expansion of Wigner-rf that only involves the d!'^^i{-K /2) symbols (see 
ADpendix lA2l for details): 



^ nm" mm' 



V(2/ + 1)(2Z" + 1) 



47r 



M" = -l'' 

where the coefficients, 



i" I' 

E/ d^m"M" (^-^ d\l"0 E/ d^mM' d\ii^i ^— ^ ^ 



^(M + M' + M")(m" -m + n)m' 



I — 'm-i m2m^ 



(-1)™! r 



im2 <^ 



sin 6d6 e 



imifl -im3p(q) 



(34) 



(35) 



have to be precomputed for a given scan strategy. Since the above expression [Eq. | |34H provides just an alternate form for the most general 
expression for bias, the discussion on different special cases given after Eq. l l28b also holds here. Thus, in case I [p(q) = 0(6) + ^{(l))] ths 
S coefficients would take the form 

r2n 



' — 7712 m3 



-1)' 



Jo 



sined6le""i''e-™^®'''\ 



(36) 







which are likely to contribute significantly for a highly constrained set of 7711,7712,7713 values and in case II (equal declination scans) 
2^177121713 contributes only for 7712 — 0, so the expression for the bias matrix becomes 



Aw = B'l 



{21' + 1) 



16tv 



I i' 

E E 

n= — l m.— — l' 



E V2l" + 1 Uv 

{ 777 — 77 ) ^ 



I" 



E ci(„_„)A/" i-^j d\j„a (- j E d^riM i-^j rfjV/O (-Jj) X 



E/ '^mW il^) E/ Pl'm' d\iim' {—j T, 



m' {M+M' + M") 



W)] 



(37) 



Here we have use the same V coefficients as defined in Eq. l l33b . In most of the cases, the beam pattern and the scan strategies have trivial 
symmetries, where only the real parts of the V coefficients contribute. 

The case of non-rotating beams, p(q) = 0, has been explicitly worked out in Appendix IbI This is a specific example where the real part 
of the r symbols can be written in a closed form and the result is: 



5R[r„,A,[p(q) =0]] - f„ 



(_^)(m'±i)/2^^2 if 771' = odd and TV = ±1 
(-l)'"'/2 2/(1 _ N^) if both m',N ^0 or even 
otherwise. 



(38) 
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So in order to compute the bias matrix for non-rotating beams, only the F symbols in Eqs. l l32t or \37i have to be replaced by the above 
/ symbols. We thus arrive at the important conclusion that the computation cost for estimating the bias for equal declination scans is the 
same as that for non-rotating beams (only the precomputed F coefficients for that scan strategy would be required). It is also reasonable to 
expect that the computation cost would be of the same order even for a broader class of scan strategies (e.g., as in case I) which can provide 
reasonable approximations to the real scan strategies. 



3 LIMITING CASES: CIRCULAR BEAM WITH FULL SKY COVERAGE & INCOMPLETE SKY COVERAGE 



The analytic expressions for bias derived in the previous section reduce to the known analytical results for non-uniform sky coverage with 
circular beams studied in Hivon et al. (2002) and our earlier results for the leading order effect of non-rotating noncircular beams with full 
sky coverage (Mitra et al. 2004). 

The special cases of circular beam and/or complete sky coverage limits are readily recovered from our general expressions. 

First, the simplest case of complete sky coverage with circular beam limit can be obtained by substituting Uim = V^Sio and /3i,n ~ 
Smo - We show in the Appendix IC 1 I that we get back the well known result 



(39) 



Hivon et al. (2002) formulated MASTER (Monte Carlo Apodized Spherical Transform Estimator) method for the estimation of CMB 
angular power spectrum from 'cut' (incomplete) sky coverage for circular beams. Substituting the circular beam limit [PimSmo] in the 
expression for bias matrix we recover the MASTER circular beam result in Appendix lC2l 



Aw = Bf 



21' + 1 
An 



i+i' 

E (2r' + i) 

'=|;-i'| 







(40) 



whereW,. = EL" = -i" \U,,'m>'\V {21" + 1). 

Finally, in Appendix IC3I we recover the general formula for leading order correction with full sky coverage for noncircular beams 
presented in Mitra et al. (2004). We substitute Uim = V^^Sio in the expression for the bias matrix and get back 

2 



Aw 



b; 



{21' + 1) 



min(l.l') 

E 

--min(l,l') 



l' 

E A'm' / dcosed'^o{e)d'^^,{e) 



(41) 



Note that due to a somewhat different definition of bias matrix in Mitra et al. (2004), for C; = [l{l + 1) /{8it^)]Ci, the results differ by a 
factor of [I' {I' + l)]/[l{l + 1)] from Eq. (38) of Mitra et al. (2004). Unfortunately, the complex form of the final expression for leading order 
correction to bias matrix with non-rotating beams presented in Eq. (43) of Mitra et al. (2004), does not allow an explicit comparison of the 
results term by term. 



4 NUMERICAL IMPLEMENTATION 

The main motivation for deriving the analytic results is to evade the computational cost and the time taken to estimate the bias using end 
to end simulations. The present work takes us one step ahead. The analysis framework we have developed can now estimate the effect of 
non-circular beams including the effect of non-uniform sky coverage for any given scan strategy. However, the numerical evaluation of 
the algebraic expression for the bias matrix derived in section [2] is also a computational challenge. As discussed in section |2] for a broad 
class of scan strategies, which often provide good approximations to real scan strategies (e.g., WMAP i23V), the computational cost for bias 
estimation can be significantly reduced to the same order as needed for equal declination scan strategies. More importantly, in situations 
where a iso-declination scan strategy only provides a crude approximation (e.g., Planck), the resulting bias estimate can be useful to judge 
the severity of the asymmetric beam effect and, hence, to decide at which level of rigor the problem needs to be addressed in the analysis 
pipeline development. 

In this section we describe the detailed computation scheme for equal declination scan strategies and show that the computational 
cost is within reach. We also suggest certain criteria for the choice of the masks in order to reduce the computational burden. Though 
the implementation of our analysis can be computationally expensive in the most general cases, even then the importance of this method 
can not be underestimated. This analysis can illuminate several "shortcuts" in the end-to-end simulations to reduce the computational cost; 
furthermore, it has the potential to eventually replace the end-to-end simulation entirely. 



4.1 Fast computation of the bias matrix 

Our scheme for fast computation of the bias is based on the alternate form of the bias expressed by Eq. l l37t . Possibility of fast computation 
of the form given by the combination of Eqs. i26i and ( I28t is being considered. 
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The final analytic form of the bias matrix for equal declination scans, given by Eq. i37b . contains infinite summations. These summations 
have to be truncated for given accuracy goals using reasonable physical insights. Let us denote the (/, m) cut-offs for the mask and beam 
by (^mask, JTimask) and (/beam, JTibcam) respectively. The choice of the numerical values for these cut-offs will be provided in the numerical 
results section. 

Computation of the final expression by naively implementing the analytic expression given by Eq. ( 137b is expensive. Three major 
innovations have been introduced in order to numerically evaluate the bias matrix: 

I. We used a smart implementation of the hierarchical summations to reduce the computational cost by a few orders of magnitude. To 
calculate three coupled loops of the form 

N N N 
i = l j = l fe = l 

naively A'^"^ operations seem necessary. However, if we calculate the summation in the following order: 

N 

V{m) := ^/(m + fc); m = 1, 2, . . . , 27V (43) 

JV JV 

S = Y.Y.V{i+j) (44) 

i = l j = l 

we effectively require just 2N'^ + A*"^ = operations. The computational gain is For A'' — 3000, this factor is 1000. This example 
is for a very simple case where all the summations have the same limits, but clearly this can be extended to the case of summations with 
unequal limits and match our analysis (See Appendix IdI for details). Here the summations within the modulus symbols in ( I37t (that is for 
each set of /, ?n, n) are computed in three stages: 

• Step I: 

V'{N) ^ J2 dmA/'d) E (I) r^'(M'+iV)[pW] (45) 

AI' — — I' — ~™bc!am 

A'' runs from —{I + m,nask) to +{l + mmask) 

• Step II: 

V\M") = Yl (f) dU (f) V\M + M") (46) 

• Step III: 

i"=0 

^mask 

E <™-n)M" (I) d';;.„ (J) v\M") (47) 

For Zbcam ~ /max the Computation time with naive implementation would scale as ~ (8/3)(2mmask + l)(2mbcam + l)/max'mask' whereas 
the above algorithm reduces the computation cost to ~ (4/3)(2mmask + l)(2mbcam + l)/max' providing a speed-up factor of ~ 21'^^-^^^. 

Mildly non-circular beams, where the BDP Pim at each / falls off rapidly with m, allows us to neglect j3im for m > mbcam- For most 
realistic beams, mbeam ~ 4 is a sufficiently good approximation (Souradeep & Ratra 2001) and this cuts off the summation over BDP in the 
bias matrix Aui . 

Soft, azimuthally apodized, masks where the coefficients Uim are small beyond m > m,nask, similarly allows us to truncate the sum 
involving Uim - Moreover, it is useful to smooth the mask in /, such the Ui,-a die off rapidly for / > /mask too. 

Clearly, small values of mbcam and mmask lead to computational speed up. Detailed discussion on mask making is shown in section l431 

II. The Wigner-d functions with argument 7r/2 occur too frequently in the above evaluation. So one possibility to reduce computation cost 
was to pre-compute all the Wigner-d coefficients dl^„ ('r/2) at once. But for / ~ 1000 this scheme is limited by disk/memory storage and/or 
program Input/Output (I/O) overhead. 

However, we may observe that, in each step of computation described in Eq. ( 147 1 only one value of / occurs in the d symbols with the 
same argument 7r/2. Hence we use an efficient recursive routine presented in Risbo (1996) that generates all the dJ„„(7r/2) at once for a 
given value of /. This allows us to compute the Wigner-d symbols efficiently and use them as constant coefficients at each step without any 
significant I/O limited operations. 

III. There are several symmetries involved in the spherical harmonic transforms and the Wigner-d symbols, which could be utilized to get 
more than an order of magnitude reduction in the computation time. 
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Figure 2. The KP2 original mask (left) and tlie azimutlially smootlied mask reconstructed from tlie KP2 mask (right). 



IV. Finally, we know that the bias matrix is not far off from diagonal, because the beams are mildly non-circular (see the plots of bias 
matrices in iSlV). So we need not compute all the elements of the bias matrix. Rather, a diagonal band (could be of triangular shape) of 
average "thickness" Al can be used to calculate the C'l estimation error with a fairly high accuracy. This will give an additional speed-up 
factor of ~ l/Al. 

While speeding up the numerical implementation of the above analysis is under progress, the estimate of computational requirement for 
the basic code for /max = 3000 is quite promising, well within the computing resources available to the CMB community (see sec. l4.3l for 
details). 

It is also interesting to note here that for different models of the same beam, which would at most differ at a highly constrained band 
in the harmonic space, the bias computation needs to be repeated only for those harmonic components. This would save large amount of 
computation. This advantage may not be available to pixel based end to end Monte Carlo simulation methods for estimating bias. 

4.2 Constructing azimuthally apodized masks 

The temperature anisotropies observed by any detector are combinations of CMB as well as foreground. The dominant contribution of fore- 
ground arises from the galactic plane. While methods of foreground removal using multi-wavelength observations exist, there is significant 
residual along the galactic plane to require masking of that region prior to cosmological power spectrum estimation. The mask is designed to 
remove the effect of regions of excessive galactic emission and spots around strong extragalactic radio sources (Bennett et al. 2003b, Saha et 
al. 2006). 

The effect of masked map on the angular power spectrum estimation has been described in the literature (Hivon et al. 2002). However, 
as we have shown, the effect of cut sky becomes nontrivial for a non-circular beam function. The sum over all the Uim modes of the mask 
is responsible for the large additional computation cost. In particular, our computational cost estimate in the previous subsection shows that 
a mask that allows us to choose a modest value of rrimask leads to a proportionally smaller computational cost. A mask whose transform is 
such that power in high m modes for a given / is suppressed would clearly serve this purpose (Souradeep et al., 2006). 

To achieve this we propose a possible method for generating an azimuthally-apodized version of a given mask0 as outlined in the 
following steps: 

I. We compute spherical harmonic coefficients Uim of the original mask U°{q). We directly suppress the power at high m, by re-scaling 

f//,„ = exp (-m * m/[a mmask^]) * Uim (48) 

corresponding to smoothing the mask along the azimuthal direction. Where a determines the extent to which power is suppressed at the 
given cut-off value of mmask- 

II. The re-scaled Uim are transformed to make an auxiliary mask U'{q). 

III. However, it is clear that mask U'{q) would allow power from contaminated regions. We should ensure that all regions where U°{q) — 
remain zero. A simple way to do that would be to multiply U'{q) with U°{q) in the pixel space. Hence, we define the final apodized mask 
as 

U^{q)^[U'{q)rU''{q), (49) 

^ For concreteness, we consider the example of Kp2 mask used in the WMAP analysis Here, for simplicity we do not consider the excised point sources. As 
this mask has also 0.6 degree radius cut around 208 locations of the point sources we first fill them up, except few, very near the galactic plane. 
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where s is a sufficiently large number (a = 1, rrimask ~ 10 and s = 12 in our example shown in figures|2]&|3) to ensure that the edges of 
the final mask are smoothed to the required level. To be more explicit, multiplying U' with Uo would introduce steps in the mask. Raising 
U' to be a sufficiently high power ensures that the amplitude of the step is reduced. 

The final mask obtained in this method is shown in the right panel of figure|2] We extract Uim from this final mask U°'{q). We show the 
Uim of the original mask KP2 and the azimuthally apodized mask Kp2" in figure[3] Clearly, the later show very rapid decrease of a mode 
with m for a given I. In the example shown here, we have significant contribution only from the first 10-20 m modes, for a given /. The 
I [/imp also dies down with I allowing us also to put a cut-off at I — Z,nask ~ 100. A reconstructed mask following this method has the 
advantage that it reduces the computation time for bias matrix by a large factor. 



4.3 Estimate of time for calculating the Bias matrix 

We have (semi)empirically estimated the CPU time required by our codes for equal declination scan strategies with apodized masks. Since 
the computation cost for equal declination scan strategies is the same as that for non-rotating beams, we use the latter for simplicity. 

We ran our codes for maximum multipoles imax of 50, 60, 65, 70, 75, 85, 90 & 100 in a 8 processor (4x AMD Opteron Dual core 
at 2.6 Ghz) node and having 16 GB of RAM. The recorded computation times are shown in the log - log plot in figure |4] (circled points). 
We fitted these points with linear functions of the form y{x) — ax + /?. The best fit is obtained for (3 — —7.6 and a — 5.5 (dotted line), 
which indicates that the observed computation time scales as l''"^, whereas algorithmically we show that it should be l^. While the theoretical 
prediction (dashed line) is also quite close to the observed values, the discrepancy is perhaps related to the simplistic implementation we have 
used in our codes and this difference should reduce in near future as the code evolves. 

Further reduction of computation is conceivable in the future: 

I. If we take into account that the bias matrix is most significant only along a narrow strip (approximately Al = 20 elements wide on 
average and not necessarily uniform - could be wedge shaped) around the diagonal of the bias matrix, the time estimate reduces by another 
power of I. The solid line in Figure|4]refers to the plot of estimated computation time against multipole in the log scale, with the slope reduced 
by 1 (i.e. 5-1=4). When /max ^ 20, of course, there is no reduction in computation cost, which is shown by the first part of the solid curve. 

II. Using the symmetries of the Wigner-d symbols, it is possible to reducing the computation cost by more than an order of magnitude. 

If all the above modifications are implemented, the bias calculation is expected to be feasible in quite reasonable time with the computing 
facilities available/dedicated to the CMB community. For example, with 1000 dual core CPUs the bias matrix for Zmax ~ 3000, mmask = 20, 
JTiboam ~ 2 and Al = 20 should be computed in 8 weeks time. Note that, this is a conservative estimate. With advanced implementation of 
the numerics having exact algorithmic scaling, the computation time can reduce by as much as one order of magnitude. 



5 DISCUSSIONS 

The inclusion of the effect of noncircular beam leads to major complications at every stage of the data analysis pipeline. The extent to which 
the non-circularity affects the step of going from the time-stream data to sky map is also very sensitive to the scan-strategy. The beam now 
has an orientation with respect to the scan path that can potentially vary along the path. This implies that the beam function is inherently time 
dependent and difficult to deconvolve. 

In our present work, we have extended our analytic approach for estimating the leading order bias due to noncircular experimental beam 
on the angular power spectrum C; of CMB anisotropy - the analytical framework now includes the effect of incomplete sky coverage and it 
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Figure 4. Circled points are the datapoints of the sample time for running the code at different multipoles. The dotted line is the best fit curve. The dashed line 
is the theoretically expected estimate for the full bias matrix. The soHd line below is the estimated computation will be required if we consider only a diagonal 
band of width Al = 20 elements in the bias matrix elements where the effect of non-circular beams is significant. 



is no more limited to only the leading order correction. It can also incorporate the case of equal declination scans without demanding any 
change in the codes or the computation cost. Though the numerical implementation does not include the effect of most general scan strategy, 
the formalism presented in this work is valid for estimating the full bias correction using a (semi)analytic perturbative method that can 
replace the computationally costly end-to-end simulation used for current CMB experiments. This work also provides an analytic framework 
to perform the simulation steps more efficiently. 

Noncircular beam effects can be modeled into the covariance functions in approaches related to maximum likelihood estimation 
(Tegmark 1997, Bond et al. 1998) and can also be included in the Harmonic ring (Challinor et al. 2002) and ring-torus estimators (Wandelt 
& Hansen 2003). However, all these methods are computationally prohibitive for high resolution maps and, at present, the computationally 
economical approach of using a Pseudo-C; estimator appears to be a viable option for extracting the power spectrum at high multipoles 
(Efstathiou 2004). The Pseudo-C; estimates have to be corrected for the systematic biases. While considerable attention has been devoted to 
the effects of incomplete/non-uniform sky coverage, no comprehensive or systematic approach is available for noncircular beam. The high 
sensitivity, 'full' (large) sky observation from space (long duration balloon) missions have alleviated the effect of incomplete sky coverage 
and other systematic effects such as the one we consider here have gained more significance. Non-uniform coverage, in particular, the galactic 
masks affect only CMB power estimation at the low multipoles. The analysis accompanying the recent second data release from WMAP uses 
the hybrid strategy (Efstathiou 2004) where the power spectrum at low multipoles is estimated using optimal Maximum Likelihood methods 
and Pseudo-Ci are used for large multipoles Hinshaw et al. 2007, Spergel et al. 2007). 

The noncircular beam is an effect that dominates at large I beyond the the inverse beam width (Mitra et al. 2004). For high resolution 
experiments, the optimal maximum likelihood methods which can account for non-circular beam functions are computationally prohibitive. 
In implementing the Pseudo-C; estimation, we have included both the non-rotating noncircular beam effect and the effect of non-uniform sky 
coverage. Our preliminary estimate shows that the computation cost for Zmax ~ 3000 is within reach. Furthermore, equal declination scan 
strategies can be trivially included in this implementation. Our work provides a convenient approach for estimating the magnitude of these 
effects in terms of the leading order deviations from a circular beam and azimuthally symmetric mask. The perturbative approach is very 
efficient. For most CMB experiments the leading few orders capture most of the effect of beam non-circularity (Souradeep & Ratra 2001). 
Our results highlight the advantage of azimuthally smoothed masks (mild deviations from azimuthal symmetry) in reducing computational 
costs. This process is more efficient as compared to the isotropic apodization of masks (Efstathiou 2006), that suffers a lot more information 
loss at the edges. The numerical implementation of our method can readily accommodate the case when pixels are revisited by the beam with 
different orientations. Evaluating the realistic bias and error-covariance for a specific CMB experiment with noncircular beams would require 
numerical evaluation of the general expressions for Aui using real scan strategy and account for inhomogeneous noise and sky coverage, the 
latter part of which has been addressed in this present work. 

It is worthwhile to note in passing that that the angular power C'l contains all the information of Gaussian CMB anisotropy only under 
the assumption of statistical isotropy. Gaussian CMB anisotropy map measured with a non-circular beam coiTesponds to an underlying 
correlation function that violates statistical isotropy. In this case, the extra information present may be measurable using, for example, the 
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bipolar power spectrum (Hajian & Souradeep 2003, Hajian et al. 2004, Hajian & Souradeep 2006, Basak et al. 2006). Even when the beam 
is circular the scanning pattern itself is expected to cause a breakdown of statistical isotropy of the measured CMB anisotropy (Hivon et al. 
2002). For a non-circular beam, this effect could be much more pronounced and, perhaps, presents an interesting avenue of future study. 
Accounting for the noncircular beam may also be crucial in the context of non-gaussianity measurements. Saha et al. (2008) developed a 
procedure to estimate the CMB power spectrum directly from the raw CMB data without the need of foreground templates. A CMB map 
made with non-circular beam but analyzed assuming a circular beam could induce higher order correlations that could be misinterpreted as 
a primordial non-gaussian signal. This is perhaps more critical for methods that first deconvolve the primordial perturbation from the CMB 
maps (Yadav & Wandelt, 2005 and Yadav & Wandelt 2008) since differences between the actual and circular beam approximations could get 
amplified and propagated through the step of deconvolution. 

In addition to temperature fluctuations, the CMB photons coming from different directions have a random, linear polarization. The 
polarization of CMB can be decomposed into E part with even parity and B part with odd parity. Besides the angular spectrum Cf^, the 
CMB polarization provides three additional spectra, Cf^, Cf"^ and CP^ which are invariant under parity transformations. The level of 
polarization of the CMB being about a tenth of the temperature fluctuation, it is only very recently that the angular power spectrum of CMB 
polarization field has been detected. The Degree Angular Scale Interferometer (DASI) has measured the CMB polarization spectrum over 
limited band of angular scales in late 2002 (Kovac et al. 2002). The DASI experiment recently published 3-year results of much refined 
measurements (Lietch et al. 2005). More recently, the BOOMERanG collaboration reported new measurements of CMB anisotropy and 
polarization spectra (Piacentini et al. 2006, MacTavish et al. 2006). The WMAP mission has also measured CMB polarization spectra (Kogut 
et al. 2003, Page et al. 2007). Connecting for the systematic effects of a non-circular beam for the polarization spectra is expected to become 
important. Extending this work to the case CMB polarization is another line of activity we plan to undertake in the near future. 

In summary, we have presented a perturbation framework to compute the effect of non-circular beam function on the estimation of 
power spectrum of CMB anisotropy taking into account the effect of a non-uniform sky coverage (e.g., galactic mask). We not only present 
the most general expression including non-uniform sky coverage as well as a noncircular beam that can be numerically evaluated but also 
provide elegant analytic results in interesting limits which are useful for gathering insights to efficiently analyze data. As CMB experiments 
strive to measure the angular power spectrum with increasing accuracy and resolution, the work provides a stepping stone to address a rather 
complicated systematic effect of noncircular beam functions. 



ACKNOWLEDGMENTS 

SM would like to thank Council of Scientific and Industrial Research (India) for supporting his research. Part of the writing of this paper was 
carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space 
Administration. We thank Kris Gorski, Jeff Jewel & Ben Wandelt for providing the reference and a code for computing Wigner-d functions. 
We thank Olivier Dore and Mike Nolta for providing us with the data files of the non-circular beam correction estimated by the WMAP team. 
We acknowledge the fruitful discussions with Francois Bouchet, Simon Prunet and Charles Lawrence. Computations were carried out at the 
HPC facility available at lUCAA. 



APPENDIX A: EVALUATION OF JnJ'mm' 

We first evaluate JUrniimm' Eq. \26\ using use the Clebsch-Gordon coefficients and sinusoidal expansion of the Wigner-d functions for 
the general case and for the case of equal declination scans. But for efficient numerical implementation in the case of non-rotating beams, we 
employ a different strategy using sinusoidal expansion of the Wigner-d functions that involve only the d}^^, (7r/2) symbols. 



Al Approach I: Using Clebsch-Gordon series and expansion of Wigner-d 

Putting Eq (O, ((FU, ((FS} & ( IeTsT i in Eq ^ we get 

j'nJ'mm' ■= / ^^^q i'in (q) i^// „// (q) DJ^^/ (q, p(q) ) (Al) 

l + l" 



/ dnqyL(„+„//)(q)Dj„,„,(q,/9(q)) 

J4n 



i+i" 



— n — m")0 (q>p(q))-D™rn' (q,p(q)) 
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(-1) 

L+l' 



i+i" 



E^LO f~<L{n + m") 



47r 



L=\l-l"\ 
-yL' (m — n — m" ) ^L'm' L' 



where 

Xmm' [P(q)] 



E^L (m — n — m ) m L r /"\] 

^H-Ti-m")l'm ^LOl'm' X(m--n-m" )m' iPlqJJ, 



df7q I>5„„/(q,p(q)). 



(A2) 



The above algebra gives us the most general expression for bias. But, as discussed in the text, special cases, which are computationally 
beneficial, often provide good approximation to the real scan strategies. For equal declination scan strategies, p(q) = p{d), 

f-2TT 



Hence, in that case, 

jll" l' r / 1 \T 



Jo Jo 



V(2; + i)(2/" + i) 

47r 



i+i" 



'^lOV'O '^InV'm" 



(A3) 



L+l' 

CLOVm' XOm' [p(^)] 

L' = |L-i' 

and the final expression for the bias matrix becomes 



647r3 



n= — lm= — l' Z"— 



I + i" L + l' I' 

C^IOV'O ^Inl" (rn-n) <^L-7ni'm A' m' C^LOi' m' XOm' 

L = |i-i"| L' = ]L~-l'] m' = -l' 

The coefficients xLm' could also be expressed in an alternative (insightful) form using Eq. dElSb as: 



\N ,1 I ,1 



Then for equal declination scan strategies one gets 



,[p{e)] ^ 27V J2 c(miv(f) dW (f) r„,,^[p(e)] 



L' 
XOn 



where [Eq. l [33t l 

Jo 

and the J coefficients become 



• njn —im'pie) 

smodoe e , 



" nm" mm' 



, ,y.+m" . V(2Z + l)(2/" + l) 

l -t) fJm" (m — n) ^ X 

i + i" i + i' 

Ey^LO ^L{n-\-m ' ) X ^ ^L' {m — n — m ) ^L' m' 
^lOV'O^lnVfrn" <-'L(_„„m")i'™ ^^LOI'm' X 

L=|i-i"| L' = |L-i'| 



d /Vm / 



(A4) 



(A5) 



(A6) 



(A7) 



(A8) 



(A9) 



A2 Approach II: Using sinusoidal expansion of Wigner-d 

We start by plugging in Eq. dElSt in Eq. Ml\ 

^nm" mm' • mm' 
J4tt 
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/■27r 

4^ io 



/ sinededLo(e) dJl'„o(e) dLn'(e) e-'™'^(^) 



47r 



i' 



M" = ~l" M' = ~l' 

j-' + + /■'"^^g»{nW-m)^ H siu 0(19 e'^" + ^''' + 6-'""' "^""K (AlO) 

io io 

This is tlie expression for the bias in its full generality. For a specified general form of /o(q), we need to precompute the integral 

Em.m^ms - (-1)'"^™^ /'^ d<j!> 6™^ T sin SdS e^^ * e'^^^"'''** . (All) 
Jo Jo 

For equal declination scans (/5(q) = 
SO the J coefficients take the form 



ju"i' _ „ c V(2Z + l)(2/" + l) 



^^™A/ (^2) '^"O ('2 ) ^ (2) '^A/"0 (2) ^ 

M=-l M" = -l" 

I' 

rfmA/' (^-^^ d\ji^, ^ — ^ r„/(A/+A/' + A/") [P(^)]- (A13) 



APPENDIX B: BIAS FOR NON-ROTATING BEAMS 

As discussed in the text, the study of the bias computation for non-rotating beams provide a framework that is computationally equivalent to 
broader class of scan strategies. Hence we treat the case of non-rotating beams with greater importance. In this case the F symbols and hence 
the final expression can be expressed in a closed form. The algebra has been worked out below. 
Substituting p{9) = in Eq. JAlOt one gets, 

jn"i' _ . ^{2l + l){2l" + l) ' , /7^^ , /ttx 

'Jnm"mm' — ^ / ^ "nA/ \ ) \ '2 / 

A/=-i 

I' 

M" = -l" M' = ~l' 

{"K 

■ n + m+m'+m" / ,%A/ + A/" + A/' / . i(M+M+M")e /Ti 1 \ 

I (—1) / smOdOe ^ . (Bl) 

The above expression is real. The proof is given below: 

• Contribution for all of M, M' , M" = 0: 

For this term the integral of the above expression is real. Therefore, if n + m + m' + m" = even this term is real (because then the factor 
^n+m+m +m real). Whcu n + m + m' + m" — odd, which means at least one of n, m + m' ,Tn" is odd, this term does not contribute, 
since d'^o{-^/2)d'o^,{n /2) = Oifm + m' = odd (follows from Eq. (6) of §4.16 of Varshalovich et al. (1988)). 

• Contribution for not all of M, M',M" = 0: 

For each set of M , M' , M" in the above summation, there exists a set —M,—AI',—M", which converts the integral of the above 
expression to its complex conjugate. Since d5„„,, (7r/2) = (-l)'"'"'d'_ w (^/2)(-l)'+'"rfm_m' (7r/2) (see Eq. (1) of §4.4 of Varshalovich 
et al. (1988)), the Wigner-d symbols give a factor of So, if n + m + m' + m" = even, the sum is real, as well as the 

factor +"1 and both are imaginary if n + m + m' + m" = odd. 



Therefore, the full summation is always real. 
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Following the discussion on the reality of the expression and using Eq. OSt we can write 



^ nm mm 



^/(2; + i)(2/" + i) 

" {m — n)m" \ ^ } r, 
I 



I" 



d'r.M [^) d\l0 i^) Yl 



M=-l 
I' 

where, as defined in Eq. l l38t . 



■(m-n)M" 



^)4/"o(^) X 



2 



du'm' [-^J fm'{M+M' + M"), 



m'N 



:= K[r„,Ar[p(q) =0]] = K 



^0 



(_;^^(m'±i)/2^y2 if m' = odd and iV = ±1 
(_l)™'/2 2/(i _ N^) ifboth7n', AT = Ooreuen 
otherwise. 



(B2) 



(B3) 



(B4) 



Note that, for "symmetric" beams = for m = odd, so in the final expression terms with m' — odd shall not contribute, that is, for 
symmetric beams fm'N contributes only when both m', A'^ = or even. 

We can now write the full expression for the bias matrix for non-rotating beams with incomplete sky coverage in a closed form as: 



Aw 



= B, 



2 {21' + 1) 



IGtt 



I i' 

E E 

n— — i m — — l' 



Y V2Z" + 1 Ui 



I" (m — n) 



(B5) 



i" 



E d[^_„^]^j„ ( — j d\iuQ ( — \ ^ djji,/ f — j d5,/o (-^ 



AI" = -l'' 
I' 



E] '^mA/' f-ir) E/ Pl'm' dM'm' (-^j .fm'{M + M'+M") 



It is obvious that the above expression is identical to the expression for equal declination scans [Eq. jilW. only the V coefficients have been 
replaced by the / coefficients. If we had started by substituting p{9) = in Eq. J28t , the above expression would take the form of Eq. l l32t 
with, again, the V coefficients replaced by the / coefficients. Which means that, the computation cost for non-rotating beams and equal 
declination scan strategies are the same, only the precomputed T coefficients for that scan strategy have to be supplied. 



APPENDIX C: CONSISTENCY CHECKS 
CI The full sky and circular beam limit 

In this appendix we recover the special case of circular beam and complete sky coverage limit. From Eq. M2\ . the full sky limit [Uim\f^5iQ} 



is obtained by replacing (7i''(m-n) with v A-KSin^8mn', and for the circular beam, we replace the BDP /3i'„j' with (5„/o. So, using the definition 
given in Eq. J38t , 



Aw = B\ 
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E 
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E '^0^ (f ) '^^o (f ) -^0 



l+l' 
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From the relation C^^^qq ~ SacSa-, (Eq. (2) in §8.5.1 of |47)), we can reduce C'qqo and C'"oo to unity, and get: 

2 



Aw = B. 



2 21' + I 



min(l,l') 



l+V 



Cl-nl'n ClOl'O rfoAT 



dj\ 
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n=-min(l,l') L'^\l~l'\ 



To get the value of X]]v=-l' ^on (f ) d'^o (f ) foN, we have to start a step back. 
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d^o ( 2 j e 



(C3) 



From the relation 



I Atv 

it follows that 



21 + 1 

1 

4tt 



E i-l^dl^ (J) (I) / e'^'sinOde f e-'-^o 



= V47r5;o5mO- 



(C4) 



The last integral J e gives 2-K5mo- So, equating out both sides and rearranging, we have 



iNe ■ , 

e sm ( 



~im<h I / 

e aq) ■- 



^/2^TT 



(C5) 



The l.h.s. is identified with the relation X]]v=-l' '^ojv (f ) t^Afo (f ) /oiv> and hence Ea. JC2b reduces to 



2 2r + 1 



min(;,i') 

E 

ri — max( — — 
min(;,i') 



l + l' 



'-■l-nl'n '^lOl'O 



L'=\l-l'\ 



V2L' + 1 



Sl'o 



00 1 2 



— -B; {21 +1) iC'f_„;/„ C;o]'0 

ri — max( , — Z' ) 

We know from Eq. (1) of §8.5.1 of Varshalovich et al. (1988) 

L-aab/Jl-J-J -J====. 

\J2a + 1 

Hence, Aiy finally reduces to the well known result 



(C6) 



Aw = bUw- 



(C7) 



C2 The circular beam limit with incomplete sky coverage 



We will show in this appendix that our formulation reduces to the analytic limit of the MASTER method of Hivon et al. (2002) for the 
incomplete sky coverage taking circular beams. Putting = V47r(5;o, as done in appendix [CT] we get: 



Au 



BfK^t E 

n— — Z 711= — I' 



16-iv 



i+i" 



E ^2FTTt/i"(,„_„) E ckw'oOi 



Lm 

Inl" (m — n) 



L + l' L' 

^ CL-ml'mCl^QliO cJoJV (^■g) '^JVO (^■2) f^'' 

L' = \L-l'\ N=-L' 



Using Eq. l lCSb . we get 



Au' = Bf'l±l± E 

n= — l m— — L' 



An 



E ^/2Fn^7^ 



l+l" 



00 



, r' 

I'm. ^ 



00 



!"=0 
2 I 



A-IT E E 



An 



n— — l m— — l' 



'■^I0l"0 '^Inl" {m-n)'^ L{-m)l' m '■^LOl'O 

2 



l+l" 



E ^2FTTf/i,.(„_„) E ^w?"oC, 

l"=i) L=\l-l"\ 



Lm 

Inl" {m — n) 



(C8) 



To arrive at Eq. A31 of Hivon et al. (2002), we first replace (m — n) by m" and then open up the modulus square. The symbol C\^ii^^ 
contributes only when m" is equal to (m — n) and also I" satisfies the triangle inequality. 

2 



Aw 



An (21' + 1) 



/ + !) E E 



n= — l rn^ — V 
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R2 ' '' 

^ V 

47r(2;' + l) ^ ^ 



^2 ^'i' ^ ^ 

5^ f^i'/m'/ Uin^n Ci„in^-,-^n Ci„i 

m'2 — — l'^ m'^ — — l'^ 



B? 



00 00 



47r(2Z' + l) ^ 



I I' 



m'^^ — l'^ m'-[ — ~l'-[ n— — l m= — l' 



(C9) 



The last summation EL-; EL=-i' C^L™/™'/ C^ini„4' simplifies to {21' + 1)/{21'I + 1)<5,,,,,,5„»^» by Eq. (5) of § 8.7.2 of Varshalovich 
etal. (1988). So, we have 



l" = \l-V\ m" = 

^!l±i (21" + 1)' ^ ' 



47r 



i"=\i-i'\ 



' I" 




(CIO) 



where W;// = EL" = -!" !t^!"m"P/(2/" + l).This matches the final expression of Hivon et al. (2002) (see Eq. [A31]). 



C3 The full sky limit with non-circular beam 



The full sky limit to the final expression should reproduce the result obtained in Mitra et al. (2004). We substitute UimV^Sio [=^> 
Ui"(^m-n) = \/47r<5i"o5m,i] in Eq. ([32) and get 



Aw 



2 {21' + 1) 



miii(l,l') 



E 

771— — min(l,l') 

L l' 

,L 
ON 



l+l' 

1000 '-^imoo y^ C*; 

L=\l-l'\ 



CIO ^Im 
1000 '-'i 



^0 , V 

l~ml'm 



^OA^ \~^) 5^ Pl'm' Ci^y^f dj^^-j^r f — j /„ 



(Cll) 



Using the relation C^^qq — SacSa'y (Eq. (2) in §8.5.1 of Varshalovich et al. (1988)) we may write Cjooo = CII^qq — 1. Then rearranging 
terms, we may write 



Aw 



B^ 



2 {21' + 1) 



lin(l.l') 

E 



777— — min(l,l') 



!' l + l' 

^ ^ A'tti' ^ ^ Ci^rnl'm Cioi'm' ^ 

m'— — V L— |Z — Z'l 
2 



E '^0^ (2') (2) 



m'N 



Using the definition of fm'N [Eq. l l38H and the expansion formula for Wigner-d [Eq. I IE15H we may write 



E '^O'^ (?) (f) -^'"'^ ^ / dcosO do^,{e). 

N=-L ■'-1 

Then, using the Clebsch-Gordon series [Eq. ( IF5H we get 



Cl-ml'm,dom'{^) Ci'i^rn' = 1 )™' ciLo (^) '^L777' (^) ■ 

Finally, putting everything together, we get the expression for the bias matrix in the full sky limit with non-circular beam as 



Aw = B, 



2 (2r + 1) 



min(l,l') I' 



E E 



m — — min(l,l') m' — — l' 

which matches Eq. (38) of Mitra et al. (2004). 



dcosO d'-^o{e)ct^^, {6) 



(C12) 



(C13) 



(C14) 



(C15) 
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APPENDIX D: FAST COMPUTATION OF BIAS MATRIX FOR NON-CIRCULAR BEAM IN CMB ANALYSIS 

Computation of the bias matrix as given by Eqs. ( |32b or l |37| ) in a naive way is too costly. For fast computation of bias using Eq. l |37| l we 
employ a smart algorithm as described in section [43] The details of the algorithm and cost estimation are given below. Possibility of fast 
computation of bias starting from Eq. \32i is being explored. 

The full expression for the bias matrix for equal declination scans, p(q) = p{d), [Eq. l|37H: 



n— — l m= — l' 
m i n ( m £^ g ]j , 1 " ) 



E 



Af" = -min(m„a3k.l") 



:"=o 

I 

d\i"o (I) E ^"Ai (I) d\m X 



(Dl) 



{m — n)M 



l' min(mi-,cam 'l') 



m ' — — m i 11 ( m ^ ,-, g^n-^ , 1 ' ) 



(2) 



The following sequence of calculation is computationally cost effective. V^''^''^ have been used as intermediate arrays. This prescription is 
only for the loops inside the modulus, so for each set of 1,1' ,m,n all the three steps have to be performed. 



• Step I: 



min(mbcam.l') 



d-rnM' (yl^) 5Z Pl'm' ^A/'m' ( — ) F, 



m' — — min(m^yaj^ ,1') 



m.'(M + M' + M") 



(D2) 



N runs from -{I + ^mask) to +(/ + Zmask)- 

• Step II 

%l [M",n, m] = J2 '^"A^ (f ) '^"0 (I) 'KJ' [M + M"; m]. 

M=-l 

• Step III 



min(m^g^gj^,l ) 



I" {m — n) 



{m~n)M" (^) c'a/"o (0 V;?/[A/",7i,m]. 



(D3) 



(D4) 



i"=0 Af" = -min(m,„j^gk.l") 

Required number of cycles to compute for each pair m, n (for Zmax ^ '"max): 

i"mi,x 

{2(/ + /"max) + 1}(2/' + l)(2m;„^ + 1) + (2/"max + 1)(2/ + 1) + X] (2^" + 1) 



(D5) 



As mentioned earlier we are interested in the total number of computation cycles in the limit /max ^ /mask. Before proceeding further we 
note that I7;"m" = is limited to only mmask modes for each /". Here mmask > 0. Then the condition for non zero 

becomes \m — n\ < mmask. This in turn implies that m — n < mmask when m — n > 0, and — m + n < mmask when m — n < 0. Then 
we see that for each n, m can run only from n — mmask to n + mmask for a total of 2mmask + 1 values so that Uinm-n are non zero. 
Thus considering two outer loops over m, n total computation cycles becomes 



^ ^ (2/ + l)(2mmask + 1) {2{l + /mask) + 1}(2/' + 1) (2m:„ax + 1) + 

1 = 2 l'=2 

^mask 

(2/Lx + l)(2/ + l)+ ^(2/" + l) 

i"=o 

(2mmask + 1) E E(2^ + 1) [{2(' + 'mask) + 1}(2/' + l)(2m:„,, + 1) + 

1 = 2 l'=2 

'mask 

(2/mask + l)(2/ + l)+ ^(2/" + l) 



(2mmask + 1) E E (2/ + 1) [2/mask(2/' + l)(2m;,,, + 1) + 
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(21 + l)(2l' + l)(2m;,,, + 1) + (2Cax + mi + 1) + E (2'" + 1) 



(D6) 



The computation cycles will be decided by the maximum power of the largest term in the above expression. Clearly the second term in the 
bracket will give the maximum contribution as it contains highest powers combined from 1,1' . Hence the total number of cycles is 



In 



(2mn,aBk + 1) E 4/" (2«')(2m;,,, + 1) = {2m^^,^ + l)(2m:,,,, + 1) x 

(2Z„,ax + 1) 



1 



^max(^max ~f" 1) 



1 



(D7) 



For iniax >• 1 the computation cost scales as (4/3)(2mniask + l)(2mJnax + l)'m. 



APPENDIX E: EXPANSION OF WIGNER-D FUNCTION 
El Motivation 

This derivation is motivated from Eq. (10) of §4.16 of Varshalovich et al. (1988). However, the motivating equation had certain inconsistency, 
as it predicts -D'„„/ {4>, 6,p) = 0ifm + m' is odd, which, in general, is not true. We rectify the formula by "reverse engineering". We start 
with the second expression of the above mentioned equation [see below for steps]: 

I 



E [dLmi i^, 0, 0) (O, J, O) Dij^jsi, (0, 9, 0) x 

M3M4 (O'f 'O) ^M4.n' (0,0, 



I 



E (0, |, 0) Dif,M3 {0: 0, 0) i^Um' (0, f , 0)] e-™'" [Step 1] 

-rm^ J2 [d^^^^ (0, |, 0) Dij^^, [e, |, 0)] e-™'" [Step 2] 



M2 = -l 

= e-'-'^DL™' e-™'" [Steps] 

= -D5„m' (|+0,7r-6», (El) 

E2 Details of the Steps in the above derivation 

• Step I: 

From Eq. (1) & (2) of §4.16, pg.ll2 of Varshalovich et al. (1988). 

D'^^,{^,0,0) = e-""^ 0^^,(0,0,0) (E2) 
dL„,{0,0,p) = Z?J„„,(0,0,0)e-'"'''' (E3) 
D'^^, (0,0,0) = (5w (E4) 

• Step II: 

From the "Addition of Rotations" formula in Eq. (3) of §4.7, pg.87 of Varshalovich et al. (1988). 

I 

J2 [DUi{c^,0i,l)Dij^,{-j,e2,p)] = Dl^,{<P,ei+62,p). (E5) 

Another way is to combine the first two remaining D symbols using Eq. (1) of §4.16, pg. 1 12 of Varshalovich et al. (1988) and then evaluate 
the following in Step III using the "Addition of Rotations" formula similar to the present method: 

E K"3 (0, f , ^) Dij^r^, (0, J, 0)] e-™'". (E6) 

M3 = -l 

• Step III: 

From Eq. (1) of §4.7, pg.87 of Varshalovich et al. (1988) we may write 

I 

J2 [D'mM{0,^,0)Di,^,ie,^,0)] = Dl^,{a,p,y) (E7) 
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where a,/3,7 are to be obtained using Eq. (66)-(70) of §1.4, pg.32 of Varsfialovicfi et al. (1988). Note that the arguments of Ihe first D 
symbol have been denoted by 02 , /32 , 72 respectively and not by ai , /3i , 71 . 

From Eq. (66) of § 1.4, pg.32 of Varshalovich et al. (1988), since 0<a<27r,0</3<7r,0<7<27r 

cos a = => a= — or ^ (E8) 
cos /3 = - cos e P = % - e (E9) 

C0S7 = 7 = (ElO) 

From Eq. (67) of § 1.4, pg.32 of Varshalovich et al. (1988) 

sina = sin7 = ^ = 1. (Ell) 
smt^ 

Combining the above equations we may write 

a = |; p ^ 7V - 9; 1 ^ \ (E12) 
E3 Final expression 

We can modify Eq. JElb by changing </)^(^— |-,&^7r — S,/9^p— ^to reach the desired expansion: 

Dl^, {(t>,e,p) = e-'"-('t'-^/2)^-^r^'(P-n/2) ^ ^gj3^ 

Yl (0, |, 0) - 0, 0) Di,^^, (0, |, 0)] . 

M2,M3 = -l 

Then using the definitions of Wigner-d functions from Eq. (1) of §4.3, pg.76 and Eq. (1) of §4.16, pg. 1 12 of Varshalovich et al. (1988), we 
get 

I 

(E14) 



This also means 



rfwC^*) = j'"^™ Y (0 e'"*'^Mm' (0] ■ (E15) 

The coefficients d}^^, (7r/2) can be directly calculated using Eq. (5) of §4.16, pg. 1 13 of Varshalovich et al. (1988) 



^-'"'(2; - 2^\/(/ + m')!(/-m')! ^^^^^ 

max{Z + m' , l — m} , i \ I 1 

(-1)'= ' > ' 



fc— max{0, m' — m} 



k j \ k + m — m' 



APPENDIX F: USEFUL FORMULAE 

• Important relations [Eq. (l)s of §4.3, §4.17 & §5.4 and Eq. (2) of §4.4 of Varshalovich et al. (1988)] 

■Dw(q,p) = e'"^''' d'^^,{9)e~"^ '' (Fl) 



yU^) = \j^^DUfl,P) = \l^!^e-^"'''dUe) (F2) 
-Dw(q,p) = (-I)'"-™ D'_„_„/(q,p) (F3) 
ii:;,(q) = (-l)'"y,_,„(q) (F4) 

Note that, unlike Mitra et al. (2004), the argument of the Wigner-d function is 9 (standard definition) not cos 9. 
• The Clebsch-Gordon series: 

Expansion of the product of two Wigner-D functions [Eq (1) of §4.6 of Varshalovich et al. (1988)]: 

-DjAini(q,/5)-Dm2n2(4/5) = 

h+h 

E„l{mi+m2) p.1 i~ \^i(ni+n2) ir;z\ 

^iimi!2m2 -'^('"1 +"12 ) ("1 +"2 ) VHi ni I2n2 ' ' 

I = |il~l2l 
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where Cj^^^j^m, the Clebsch-Gordon coefficients. 

The special case of spherical harmonics [Eq. (9) of §5.6 of Varshalovich et al. (1988)]: 

Yi,^A^)Yi,m,{fi) = (F6) 



h+h 

E 

i = |il-i2| 



El {2ll + l)(2/2 + 1) ^iO ^i(™i+m2)v 
y 47r(2Z + 1) ^hOhO ^Iimil2m2 +"12 ) I^IJ ■ 



In modifying the above equations (from Varshalovich et al. (1988)) we have used the fact that the Clebsch-Gordon coefficients Cj'^^i^^^ 
vanish if m 7^ mi + m^. 
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